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1 Introduction 

Convection of an electrically conducting fluid in a rotating system represents a basic 
dynamical process in planetary interiors and in stars. Astrophysicists and geophysi- 
cists have long been interested in the mechanisms that govern the convective heat 
transport and the generation of magnetic fields by convection in those systems. The 
availability in recent years of large scale computer capacities has permitted numeri- 
cal simulations of detailed models for those processes. Only the larger length scales 
can be taken into account in those computations, of course, and eddy diffusivities 
are usually introduced to model the influence of the smaller unresolved scales of the 
turbulent flows. 

A difficulty arises from the fact that it can not generally be assumed that all eddy 
diffusivities are equal. First, they apply to scalar as well as to vector quantities, 
such as temperature and magnetic fields. Secondly, the diffusivities in the absence 
of turbulence differ enormously such that the turbulence may not be sufficiently 
strong to equalize them. In the Earth's liquid core, for instance, the magnetic diffu- 
sivity is large enough to be taken into account without the consideration of an eddy 
contribution while a comparable eddy viscosity would have to exceed the probable 
molecular viscosity value by a factor of at least 10 6 . As has been demonstrated in 
the past ifTTI the dynamics of convection in rotating spherical shells and its dynamo 
action are very sensitive to ratios of diffusivities, especially to the Prandtl number 
around its usually assumed value of unity. 

More complex methods for treating effects of turbulence could eventually be 
used, such as k — e-models, and the undoubtedly important anisotropy of turbulent 
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diffusivities in rotating systems could also be considered. Since these influences are 
difficult to evaluate, however, we shall restrict the attention in this paper to a min- 
imum of physical parameters. Inspite of this restriction a number of characteristic 
features of convection and its dynamos in rotating spherical shells can be demon- 
strated that are likely to exist as coherent spatio-temporal structures in natural sys- 
tems. 



2 Mathematical formulation of the problem and methods of 
solution 

We consider a rotating spherical fluid shell of thickness d and assume that a static 
state exists with the temperature distribution T$ = Tq — [5d 2 r 2 /2. Here rd is the 
length of the position vector with respect to the center of the sphere. The grav- 
ity field is given by g = —dyr. In addition to d, the time d 2 /v, the temperature 
V 2 /yad 4 and the magnetic flux density v(jlip) 1 / 2 jd are used as scales for the di- 
mensionless description of the problem where v denotes the kinematic viscosity of 
the fluid, fc its thermal diffusivity, p its density and jj. is its magnetic permeability. 
Since we shall assume the Boussinesq approximation material properties are re- 
garded as constants except for the temperature dependence of the density described 
by a = — (dp / dT) / p which is taken into account only in the gravity term. Both, 
the velocity field v and the magnetic flux density B, are solenoidal vector fields for 
which the general representation 

v = Vx (Vaxr)+Vwxr, B = V x (V/i x r) + Vg x r, (1) 

can be employed. By multiplying the (curl) 2 and the curl of the Navier-Stokes equa- 
tions of motion by r we obtain two equations for u and w, 

[(V 2 -d r )j^2 + Td,p]V 2 M + T^w-^20 = -r-Vx [Vx (v-Vv-B-VB)], (2) 
[(V 2 - 3)jSf 2 + xd 9 ]w - z£u = r • V x (v • Vv - B • VB), (3) 

where d t denotes the partial derivative with respect to time t and where <9<p is the 
partial derivative with respect to the angle <p of a spherical system of coordinates 
r,9,(p. For further details we refer to ifTTI . The operators _S?2 and =2 are defined by 

J? 2 = -r 2 V 2 + d r (r 2 d r ), 
J2 = r cos V 2 - (jSf 2 + rd r ) (cos Gd r - r" 1 sin 9d e ). 

The heat equation for the dimensionless deviation from the static temperature 
distribution can be written in the form 



V 2 @ = P{d, + v • V)0, 



(4) 
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Fig. 1 Convection columns in a rotating spherical fluid shell for x = 10 4 ,fi = 2.8- 10 5 ,P = l.Blue 
and red surfaces correspond to a constant positive and negative value of the radial velocity. 



and the equations for h and g are obtained through the multiplication of the equation 
of induction and of its curl by r, 

V 2 jSf 2 A = P m [d,^2h - r • V x (v x B)] , (5) 
V 2 J? 2 g = P m [d,^ 2 g - r • V x (V x (v x B))]. (6) 

The Rayleigh number R, the Coriolis number T, the Prandtl number P and the mag- 
netic Prandtl number P m are defined by 

VK V K A 

where X is the magnetic diffusivity. For the static temperature distribution we have 
chosen the case of a homogeneously heated sphere. This state is traditionally used 
for the analysis of convection in self-gravitating spheres and offers the numerical ad- 
vantage that for Rayleigh numbers close to the critical value R c the strength of con- 
vection does not differ much near the inner and outer boundaries. As the Rayleigh 
number increases beyond R c heat enters increasingly at the inner boundary and is 
delivered by convection to the outer boundary. When R reaches a high multiple of 
R c the heat generated internally in the fluid becomes negligible in comparison to the 
heat transported by convection from the inner to the outer boundary. 
Fixed temperatures and stress-free boundaries, 

u = <9 2 u = d r {w/r) =0=0 at r = r x = tj/(1 -tj) and r = r = 1/(1 - 77), (8) 

will be assumed where 77 denotes radius ratio, 77 = r,/r . In the following only 
77 = 0.4 will be used. For the magnetic field electrically insulating boundaries are 
assumed such that the poloidal function h must be matched to the function h which 
describes the potential fields outside the fluid shell, 
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Fig. 2 Time series of energy densities of convection for P = 1, T = 10 4 and R = 2.8 x 10 5 , 
3.0 x_10 5 , 3.5 x 10 5 , 7 x 10 s , 12 x 10 5 , (from top to bottom). Solid and dashed lines indicate E, 
and £",, respectively. The Nusselt number Nit; is indicated by dotted lines and measured at the right 
ordinate. 



g = h-h = d r (h — h) =0 at r = r, = 17/(1 — rj) and r = r = 1/(1 -tj). (9) 

But computations for the case of an inner boundary with no-slip conditions and 
an electrical conductivity equal to that of the fluid have also been done 1111 . The 
numerical integration of the equations together with the boundary conditions pro- 
ceeds with the pseudo-spectral method as described in (13\ and 11141 which is based 
on an expansion of all dependent variables in spherical harmonics for the 9 , <p- 
dependences, i.e. 

" = E U ? ( r > r ) p r ( cos ) exp-j)>w<p } ( 1 0) 

l.m 

and analogous expressions for the other variables, w,©,h and g. P" 1 denotes the as- 
sociated Legendre functions. For the r-dependence expansions in Chebychev poly- 
nomials are used. For further details see also [ 1 1 1. For the non-magnetic convection 
calculations to be reported in the following a minimum of 33 collocation points in 
the radial direction and spherical harmonics up to the order 64 have been used. The 
resolution has been increased to 41 or 55 collocation points and spherical harmonics 
up to the order 96 or 128 in the case of dynamo simulations. 
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Fig. 3 Localized convection for R = 7 x 10 5 , T = 1.5 x 10 4 ,P = 0.5 The streamlines, rdu/dq> = 
const, (upper row) and the isotherms, = const, (lower row), are shown in the equatorial plane 
for equidistant times (from left to right) with At = 0.03. 



3 Convection in rotating spherical shells 

For an introduction to the problem of convection in spherical shells we refer to the 
review and to the respective chapter in the book Q ■ Additional information can 
be found in the papers by Grote and Busse @, Jones et al. Q, Christensen J4) 
and Simitev and Busse iflOl . Typically the onset of convection occurs in the form 
of progradely propagating thermal Rossby waves as illustrated in figure 1 . Only for 
low Prandtl numbers P, i.e. P < 10/ ' yfx according to UJ, the onset occurs in the 
form of inertial waves attached to the outer equatorial boundary of the fluid shell as 
is discussed in lfl31l.lfT6l.lfT1. 

Because of the symmetry of the velocity field with respect to the equatorial plane 
it is sufficient to plot streamlines in this plane, given by rdu/dr = const., to charac- 
terize the convection flow. This has been done in figures 3 and 5. Even in the case of 
turbulent convection the part of the velocity field that is antisymmetric with respect 
to the equatorial plane is rather small as long as the parameter T is sufficiently large. 

As the Rayleigh number/? grows beyond its critical value R c , the thermal Rossby 
waves become modified by a sequence of bifurcations similar to those found in 
other problems of convection. First, oscillations of the amplitude are observed, then 
another bifurcation causes a low wavenumber modulation as function of the azimuth 
ifTol . Finally, a chaotic state of convection is obtained. 
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4 Chaotic Convection 

The sequence of transitions can also visualized through the time dependence of 
average quantities such as the contributions to the kinetic energy density. These are 
defined by 



where the angular brackets indicate the average over the fluid shell and where u 
refers to the azimuthally averaged component of u and u is given by u = u — u. At 
low supercritical Rayleigh numbers the energy densities corresponding to steadily 
drifting thermal Rossby waves are constant in time and are not included in figure 2. 
The onset of vacillations manifests itself in the sinusoidal oscillations of the kinetic 
energies as shown in the top plot of figure 2. Also plotted in figure 2 is the Nusselt 
number Nuj measuring the efficiency of the convective heat transport at the inner 
boundary, 



where the double bar indicates the average over the spherical surface. E t describes 
the energy density of the differential rotation which increases strongly with increas- 
ing R as can be noticed in the lower plots of figure 2. This increase is caused by the 
strong azimuthal Reynolds stress exerted by the convection eddies resulting from 
their inclination with respect to radial direction as is apparent in figure 1 . The in- 
creasing shear of the differential rotation tends to inhibit convection, however, in 
that it shears off the convection eddies. This is a consequence of the nearly two- 
dimensional nature of the dynamics in a rotating system: Because of the require- 
ment that the structure of convection approaches closely the Taylor-Proudman con- 
dition, there is no possibility for a reorientation of the convection rolls as happens 
in non-rotating systems. In the rotating sphere convection thus generates the agent 
that tends to destroy it. A precarious balance in the form localized convection is 
the result. As shown in figure 3 convection occurs only in a restricted azimuthal 
section of the spherical shell where its amplitude is strong enough to overcome the 
inhibiting influence of the shear. For the geostrophic zonal flow it does not matter 
whether it is driven locally or more uniformly around the azimuth. Since in the non- 
convecting region thermal buoyancy accumulates, the advection of this buoyancy 
by the differential rotation strengthens and stabilizes the localized convection. 

At even higher Rayleigh numbers this balance no longer works and instead of a 
localization in space the localization of convection in time is initiated as is shown in 
figures 4 and 5. Here convection exist only for a short period while the differential 
rotation is sufficiently weak. As the amplitude of convection grows, the differen- 
tial rotation grows even more strongly since the Reynolds stress increases with the 
square of the amplitude. Soon the shearing action becomes strong enough to cut off 



£ P = (|Vx (Vu x r) | 2 )/2, E, = (| Vw x r | 2 )/2 
E p = (| V x (Vh x r) | 2 )/2, E, = (\Vwxr | 2 )/2 



(lib) 



(Ha) 




(12) 
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Fig. 4 Relaxation oscillations of chaotic convection in the case T = 10 4 ,fi = 6.5 • 10 s , P = 0.5. 
The energy densities E, (solid line), E t (dashed line), E p (dotted line) and the Nusselt number 
(dash-dotted line, right ordinate) are shown as function of time. 



convection. Now a viscous diffusion time must pass before the differential rotation 
has decayed sufficiently such that convection may start growing again. It is remark- 
able to see how the chaotic system exhibits its nearly periodic relaxation oscillations 
as shown in figure 4. 

In figure 5 a sequence of plots is shown at four instances around the time of a con- 
vection peak. At first there is hardly any convection, - the dotted lines just indicate 
zero. At the next instance the differential rotation as shown by the upper row has 
decayed sufficiently such that convection columns can grow reaching nearly their 
maximum amplitude in the third plot. At the same time the differential rotation has 
grown as well and begins to exert its inhibiting effect such that convection decays at 
the fourth instance of the sequence, while the differential rotation reaches its max- 
imum. It should be mentioned that localized convection and relaxation oscillations 
occur at moderate Prandtl numbers of the order unity or less. At higher values of P 
Reynolds stresses are no longer sufficiently powerful to generate a strong differen- 
tial rotation. Instead variations of the temperature field caused by the dependence of 
the convective heat transport on latitude induce a differential rotation in the form of 
a thermal wind. 

The convective heat transport in the case of localized convection as well as in 
the case of relaxation oscillations is much reduced, of course, relative to a case 
without strong differential rotation. This causes the magnetic field to enter the prob- 
lem in a crucial way provided the electrical conductivity is sufficiently high. By 
putting brakes on the differential rotation through its Lorentz force the magnetic 
field permits a much higher heat transport than would be possible in an electrically- 
insulating fluid. This is the basic reason that the Earth's core as well as other planets 
with convecting cores and rotating stars exhibit magnetic fields. A demonstration of 
this effect is seen in figure 6 where by chance the convection driven dynamo was 
just marginal such that it could not recover after a downward fluctuation of the mag- 
netic field. Hence the relaxation oscillations with their much reduced average heat 
flux take over from the dynamo state. 
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Fig. 5 Sequence of plots starting at t = 2.31143 and equidistant in time {At = 0.01) for the same 
case as in Fig. 8. Lines of constant Uy and streamlines rsin 8dgh =const. in the meridional plane, 
are shown in the left and right halves, respectively, of the upper row. The lower row shows corre- 
sponding streamlines, rdu/d(p = const., in the equatorial plane. 



E 




Fig. 6 Transition from a dynamo state to a state of chaotic relaxation oscillations for T = 1.5 x 
10 4 ,R = 1.2 x 10 6 ,P = P„, = 0.5. The energy density E, (dotted line), the total magnetic energy 
density (dashed line) and the Nusselt number Niij (solid line, right ordinate) are shown as function 
of time. 



5 Distinct turbulent dynamos at identical parameter values 

Convection driven dynamos in rotating spherical fluid shells are often subcritical as 
is apparent in figure 6, for instance. At somewhat higher Rayleigh numbers con- 
vection with a strong magnetic field will persist. On the other hand the dynamo will 
decay when the magnetic field is artificially reduced to, say, a quarter of its averaged 
energy. There thus exists the possibility of a convection driven dynamo state and of a 
non-magnetic convection state at identical values of the external parameters R, x,P. 
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Fig. 7 Two distinct dynamos at identical parameter values, T = 5x 10 3 ,R = 5 x 10 5 ,P = P m = 1. 
Upper (lower) plots show time series of quadrupolar (dipolar) magnetic energy densities. Thick 
lines indicate mean torodial (solid lines) and poloidal (dashed lines) energy components. Thin 
lines indicate the same for the fluctuating components. 



The bistable coexistence of a non-magnetic convection state and a dynamo state is 
typical for subcritical bifurcations as in the analogous coexistence of laminar and 
turbulent states in shear flows. 

More surprising is the fact that two different turbulent dynamo states can exist at 
identical values of the external parameters which now should include the magnetic 
Prandtl number P m . An example is shown in figure 7 where a spherical dynamo of 
mixed parity, i.e. with dipolar and quadrupolar components, and a purely quadrupo- 
lar dynamo evolve in time at identical external parameters. In the latter case the 
exponential decay of dipolar disturbances is clearly demonstrated. In figure 7 mag- 
netic energy densities have been plotted that are defined in analogy to expressions 
(12), 

M p = {\Vx{Vhxr)\ 2 }/2, M, = (\Vg t xr\ 2 )/2 (13a) 
M p = (|Vx(V^xr)| 2 )/2, M t = (| V^xr| 2 )/2. (13b) 

An even more surprising case is that of two convection driven dynamos without 
any distinction in symmetry, just with differences in the magnitude of various energy 
densities as is apparent from figure 8. While a strong mean poloidal magnetic field as 
shown in the left half of figure 8 acts as an efficient brake on the differential rotation 
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Fig. 8 Time series of two different chaotic attractors are shown - a MD (left column (a,b)) and a 
FD dynamo (right column (c,d)) both in the case R = 3.5 x 10 6 , t = 4 x 10 4 , P = 0.5 and P m = 1. 
The top two panels (a,c) show magnetic energy densities, and the bottom two panels (b,d) show 
kinetic energy densities in the presence of the magnetic field. The components M p are shown by 
thick solid black lines, while X,,X P , and X, are shown in red green and blue respectively. X stands 
for either M or E. 
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Fig. 9 The upper row shows the hysteresis effect in the ratio of magnetic to kinetic energy, E/M, 
at T = 3 x 10 4 (a) as a function of the Prandtl number in the case of R = 3.5 x 10 6 , P/P„, =0.5; 
(b) as a function of the ratio P/P„, in the case of R = 3.5 x 10 6 , P = 0.75 and (c) as a function 
of the Rayleigh number in the case P = 0.75, P m = 1.5. Full and empty symbols indicate FD 
and MD dynamos, respectively and circles and squares indicate the two hysteresis branches. The 
critical value of R for the onset of thermal convection for the cases shown in (c) is R c = 659145. A 
transition from FD to MD dynamos as P/P„, decreases in (b) is expected, but is not indicated owing 
to lack of data. The lower row shows the value Nui of the Nusselt number at r = r, for the same 
dynamo cases. Values for non-magnetic convection are indicated by triangles for comparison. 
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Fig. 10 The upper part shows the ratio between the energy densities of the fluctuating and the 
mean poloidal magnetic field, M„/M p , as a function of the Prandtl number at T = 3 x 10 4 (denoted 
by 30 in plot)and at t = 4 x 10 4 (denoted by 40 in plot) by symbols analogous to those in fig 10. 
In all cases R = 3.5 x 10 6 and P m = 2P. The lower part shows the value Nui of the Nusselt number 
for the same dynamo cases. 



as measured by E t , it also inhibits convection. The alternative dynamo on the right 
side of the figure is characterized by a relatively weak mean magnetic field and 
dominant fluctuating components. Here the kinetic energy densities of convection 
are larger, but the differential rotation is still much weaker than it would be in the 
non-magnetic case. 

The bistability in the form of two different types of dynamos is not a singular 
phenomenon, but exists over an extended region of the parameter space. Since the 
region includes values of the Prandtl number of the order unity which has been 
the preferred value of P in many simulations of convection driven dynamos the 
phenomenon of bistable dynamos is of considerable importance. The first results 
of our extensive numerical simulations can be found in lfl2l . A typical diagram is 
shown in figure 9. The extended regime of coexistence of the two types of dynamo 
is bounded by transitions where one of the two dynamos ceases to be stable and 
evolves into the other one within a magnetic diffusion time. A basic reason for the 
competitiveness of both dynamos is that they exhibit essentially the same convective 
heat transport as measured by the Nusselt number Niij. The lower half of figure 9 
not only demonstrates the surprising coincidence of the heat transports of the two 
dynamo types, but also indicates that these heat transports by far exceed those found 
in the absence of a magnetic field. 

Recently the computations have been extended to higher rotation rates as shown 
in figure 10. At z = 4 x 10 4 the same phenomenon of bistability has been found as 
at z = 3 x 10 4 . Again, the two kinds of dynamo exhibit the same heat transport as 
indicated in the lower part of figure 10. 
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The existence of two distinct turbulent states is a rare phenomenon, although exam- 
ples exist in non-magnetic hydrodynamics, see, for instance, 13 . In magneto- 
hydrodynamics an electrically conducting fluid in the presence of a magnetic field 
offers new degrees of freedom which allow more than a single balance between the 
various forces operating in turbulent states. Initial conditions thus determine which 
of the competing states is actually realized. 

The possibility of bistability could be of interest for the interpretation of plane- 
tary and stellar magnetism. Magnetic hysteresis effects associated with stellar oscil- 
lations may eventually be explained in this way. Anyone involved with numerical 
simulations of convection driven dynamos should be aware that his solutions could 
change drastically after different initial conditions have been introduced. 
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